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Abstract. Here we present complex resonance states (or Siegert states), that describe 
the tunneling decay of a trapped quantum particle, from an intuitive point of view 
which naturally leads to the easily applicable Siegert approximation method that can 
be used for analytical and numerical calculations of complex resonances of both the 
linear and nonlinear Schrodinger equation. Our approach thus complements other 
treatments of the subject that mostly focus on methods based on continuation in the 
complex plane or on semiclassical approximations. 
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1. Introduction 

The tunneling decay of particles trapped in an external potential V{yi) is a problem 
often encountered in nuclear, atomic and molecular physics. The most famous example 
is arguably the alpha decay of nuclei for which Gamow calculated the decay rates 
semiclassically [1]. 

For a particle with mass m such a decay can be described by resonance eigenstates 
ips of the Schrodinger equation 

(-|^V^ + V(x)) = £^s (1) 

with complex eigenenergy S = Es — ir/2, often refered to as Siegert resonances, where 
the wavefunction satisfies outgoing wave (or Siegert) boundary conditions [2], i. e. the 
wavefunction is given by an outgoing plane wave for |x| — )■ ±oo (cf. equation ([S])). The 
imaginary part is also refered to as the decay rate since it leads to an exponential decay 
of the wavefunction, 

^5(x, t) = exp{iSt/h)iJs{^, 0) = exp{-Tt/{2h)) exp{iEst/h)^s{^, 0) . (2) 

On the other hand, the tunneling decay of trapped particles is closely related to 
the problem of scattering of particles off the same potential since quantities like the 
transmission coefficient \T{E)\'^ or the scattering cross section show characteristic peaks 
near the resonance energies S, which can be described by a Lorentz or Breit-Wigner 
profile 

^^^^^^'= {E-Esyl{T/2y 

the width of which is given by the decay rate F [2] . 

As mentioned above, the decay rates can be calculated in the manner of Gamow 
using semiclassical approximations. Though these approximations are straightforward 
and easy to use they are often not very precise, only providing an order of magnitude. 
On the other hand there are powerful complex- scaling based methods (see e. g. Oil]) 
(including related methods such as complex absorbing potentials) where the spatial 
coordinate is rotated to the complex plane, x — xexp(i^) by some sufficiently large 
angle 6 to make the resonance wavefunctions square integrable, enabling the use of the 
usual techniques for calculating ordinary bound states. These methods are precise and 
highly efficient yet quite sophisticated and, apart from rare exceptions, only suited for 
numerical calculations. Alternative techniques, like, e.g., the stabilization method [5] 
have both some advantages and disadvantages compared to complex scaling and are, 
however, not very intuitive and only aim at numerical applications. While there are 
a number of excellent texts (see, e. g. [3llll[6] and references therein) that discuss the 
mathematical aspects of the problem (e. g. analytical continuation of the wavefunction in 
the complex plane) and the computational methods mentioned above, a complementary 
treatment assuming a different point of view could prove valuable. 
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In this article we want to draw attention to the scarcely known Siegert 
approximation method for calculating complex resonances which is both intuitive 
and easy to implement but does not rely on semiclassical arguments. It yields good 
results for narrow resonances (i.e. r/2 <^ Es-, where the resonance energy Es is 
measured relative to the potential energy at |x| — > oo) and it can in some cases 
lead to closed form analytical approximations for the decay coefficient. Another 
advantage of this method is its straightforward applicability to resonances of the 
nonlinear Schrodinger equation that occur, e. g., in the context of trapped Bose-Einstein 
condensates [TfS] since it does not require properties like the linearity or analyticity of the 
differential equation. While direct complex scaling [9|[T0] has succesfuUy been applied to 
nonlinear Siegert resonances it is much less efficient than in the linear case and requires 
substantial modifications. Complex absorbing potential methods, which have proven 
more efficient in this context [mill], also require considerable modifications compared 
to the linear case. 

In the derivation of the Siegert approximation method given here, complex 
resonances are reviewed from an alternative, rather intuitive point of view which 
complements the usual more mathematical treatments of the problem by emphasizing 
some important aspects like the role of the continuity equation in this context as 
well as the similarities and differences between complex resonance states and so-called 
transmission resonances, i. e. scattering states corresponding to the maxima of the 
transmission coefficient (or scattering cross section respectively) for real eigenenergies. 

This paper is organized as follows: In section [2]the Siegert approximation method is 
described, focussing on resonances of the one- dimensional linear Schrodinger equation 
for the sake of simplicity. In section [3] the method is illustrated by means of several 
analytical and numerical applications. Finally, the main aspects of the article are 



summarized in section HI Appendix A contains a MATLAB code for calculating 
resonances. 

2. The Siegert approximation method 

For narrow resonances (i.e. r/2 <^ E), one can often neglect the decay rate F in ([1]) 
and thus obtain an approximate resonance wavefunction V^approx and an approximate 
real part -Eapprox of the resonance energy with very little effort. 

For the sake of simplicity we first consider symmetric one-dimensional finite range 
potentials 



v{x) = i , r (4) 



1 ' 


\x\ 




\x\ 



with finite range a and V{—x) = V{x). An example of such a potential is the double 
barrier shown in the right panel of figured] (bold blue curve) which can safely be assumed 
to be approximately equal to zero for > 20. Also shown is the square |'?/^(x)p of the 
most stable resonance wavefunction which is strongly localized between the potential 
maxima and inherits the symmetry of the double barrier potential. 
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Apart from the Siegert resonances ips with complex energies Es — ir/2 obtained 
for outgoing wave boundary conditions 



ij's{±a) = ±iksi^s{±a), ks = pm{Es - ir/2)/^ (5) 

(the prime denotes a derivative by x) these potentials also possesses so-called 
transmission (or unit) resonances iPt{x) corresponding to real energies Et for which the 
potential is completely transparent, i.e. for the corresponding transmission coefficient 
we have \T{Et)\'^ = 1. We will see further below, that for narrow resonances iprix) and 
Et are good approximations to ips and Es respectively. 

To obtain the transmission coefficient |Tp we consider transmission through the 
potential V{x) which is characterized by the following boundary conditions for the 
scattering wavefunction ip{x): On the left hand side we have a superposition of an 
incoming and a reflected plane wave 

ip{x) = Aexp{ikx) + Bexp{—ikx), x < —a (6) 



where k = y2/rnE/h is the wavenumber corresponding to an energy E of the incoming 
wave. On the right hand side we only have an outgoing wave, 

ip{x) = C exp(i/i;x), x > a . (7) 

Thus the transmission coefficient reads 

\T\^ = \C/A\^. (8) 

We immediately see that for the transmission resonances with |Tp = 1 we have |C| = |y4| 
(and -B = 0), so that we can always achieve C = Ahy multiplying the wavefunction with 
a constant phase factor. Thus the boundary conditions for the transmission resonces 
il)T can be written as 

ipT = Aexp{ikTx), \x\ > a (9) 

or equivalently as 



?/;^(±a) = ikTi^T{±a) with kr = \J2mET/h. (10) 

The symmetry of both the Siegert resonance wavefunction \ips{~x)\^ = \tps{x)\'^ and the 
transmission resonance wavefunction \iIjt{—x)\'^ = |-?/'r(2^)P imply that the derivatives 
of l^psl"^ and liprl"^ must vanish at x = 0. Thus the respective boundary conditions ([5]) 
and (llUp can be recast in the form 

(IV^a(x)n'L=o = 0' tpM = ik^Ma), ae{S,T} (11) 

with k^ = y^2^^/h, £s = Es- ir/2, £t = Et- 

Therefore it is quite intuitive that for r/2 ^ Es we can make the approximations 
Es ~ Et and 

, , f i^rix) <x <a 
i'six) ~ < , / \ • (12) 

This was shown rigorously by Siegert [2]. Note that this approximate 
correspondance between complex energy Siegert resonances describing decay and real 
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energy transmission resonances, which generally holds for narrow resonance widths, is 
in fact one of the main reasons for considering complex resonances in the context of 
scattering [2]. Now we have approximations for the wavefunction and the real part of 
the eigenenergy but what about the imaginary part, i.e. the decay rate ? 

To this end let us assume that the potential V{x) has local maxima at a; = ±6 with 
< 6 < a (as, for example, the potential shown in figure ([1]) ). Then the probability 
of finding the particle described by our Schrodinger equation ([T]) 'inside' the potential 
well, i.e. in the region —h<x<+h between the potential maxima is given by the norm 

N = t |V^5(x)Pdx ^ 2 r |V^T(x)Pdx (13) 

J-b Jo 

of the wavefunction inside the well. The exponential decay behaviour of the resonance 
wavefunction given in equation ([2]) implies that the norm decays according to 

dtN = -^N (14) 
so that the decay rate can be written as 

r = -RM. (15) 

AT ^ ^ 

The time derivative dtN = j'^f^dt\ijjs{x,t)\'^dx ~ 2 dt\i/jT{x,t)\'^dx of the norm can be 
found by means of the continuity equation for the resonance wavefunction ips which 
reads 

dtps = -j's (16) 
with ps{x,t) = \ips{x,t)\'^ and the probability current 

Js = irsi^'s - ^s^s) ■ (17) 



Integrating the continuity equation f[T6l) from x = —b to x = b v/e obtain dtN = 
~{js{b) — jsi~b))- Equation f[T^ implies that the currents are approximately given 
by js{b) ~ jrib) and jsi—b) ~ —jrib). Furthermore, the current jr corresponding to 
the transmission resonance ipx does not depend on the position x and in particular 
jrib) = jxio)- For x > a the transmission resonance wavefunction is given by a 
plane wave iprix) = (a) exp (1/^^(3^ — a) so that the current at x = a simply reads 
jrio,) = hkT\ip{a)\'^ /m. Thus the decay coefficient f[T^ becomes 

JtlM^Wdx m j'^\ijT{xWdx- ' 
We call ( 1T8|) the Siegert formula. Note that it only depends on Et and iI)t so that the 
exact values Es and ifjs are not required. 

In general, the Siegert approximation method consists of two steps: 

(i) Neglect at first the imaginary part T /2 (also called decay coefficient) of the 
complex resonance energy in the Schrodinger equation for the Siegert resonance 
wave function in order to obtain approximations to both the wave function and 
the real part of the resonance energy. (In the symmetric barrier case discussed 
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above this is done by calculating the energies and scattering states for which 
the transmission probability of the corresponding scattering problem has a local 
maximum or, equivalently, directly use the symmetry of the problem expressed in 
the boundary conditions given in equation (ITT]) .) 

(ii) Use the continuity equation to obtain a Siegert formula, i. e. an approximate 
expression for the decay coefficient F (as given in equation (flSj) for the symmetric 
barrier case) which only depends on the approximate energy and approximate 
wavefunction calculated in step 

The Siegert approximation method yields good results for narrow resonances, 
i. e. whenever the imaginary part r/2 of the resonance energy is small compared to 
the real part E. Unlike common discussions of complex resonance states the above 
derivation is based on the intuitive picture of a matter wave flowing out of a potential 
well, emphasizing the role of the conservation of the probability current expressed by 
the continuity equation. 

The above treatment for symmetric potentials can be straightforwardly generalized 
to resonances of asymmetric barrier potentials where the matter wave is localized 
between two maxima at a; = 6_ < and a; = 6+ > with two flnite ranges a_ , h_ and 
a_|_,6_|_ but then the actual calculations in step (jl]) are generally more difficult because 
the symmetry condition (ITT]) no longer applies. Now one has to consider the problem 
of transmission through the barrier and calculate the states 'j/'approx and the respective 
real energies -Eapprox which correspond to local maximima of the transmission coefficient 
with |T(i?approx)P < 1- In analogy to the symmetric case one finds that the relation 
(fT8|) for the decay coefficient is generalized to 

hi fcapprox|V^approx('^+) I ~l~ h /c_a|'?/'approx('3'— ) | ( ^ (W 

'^/fe„^ |V^approx(a;)pda; 



where kapprox = Y2mi?approx/^- An important special case of an asymmetric barrier is 
a trap which is open on one side only whereas there is an impenetrable barrier on the 
other side. If the impenetrable barrier is on the left hand side (at x = a_) we obtain 
the boundary condition 

|^approx(«-)|' = (20) 

and the wavefunction ^/'approx and the corresponding energy -Eapprox can be obtained by 
finding an approximate solution to the system of equations given by ( 120|) and 

V'lpprox('^+) = i^approx^/'lpprox(^+) • (21) 



An example of such a calculation is given in section 13.21 Note that for such a potential 
exact solutions for real energies do not exist. We further note that potentials with infinite 
range can usually be approximated by finite range potentials by choosing appropriate 
values for a± and that the Siegert approximation method can be generalized to two and 
three dimensions (see [7] for a detailed discussion). 
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3.1. The finite square well potential 




Figure 1. Barrier potentials (bold blue) and approximate squared wavefunctions 
(red) of the corresponding most stable resonance states (units with h = m = 1). Left: 
Finite square well with = -5, L 3. Middle: Delta Shell potential ^ with 
A = 10, L = 1 (The Dirac delta function is symbolically represented by a vertical line). 
Right: Double barrier potential with VJ) = 1, A = 0.1. 



As a first analytically solvable example we consider the finite square well potential 

Vi^) = I 1"^! ^ ^ (22) 

^ ' \ Vo<0 \x\<L ^ ' 

with width 2L which is shown in the left panel of figure [T] As discussed at the beginning 
of the preceding section, the complex resonances for such a symmetric potential can 
be calculated approximately by finding its transmission resonances which satisfy the 
symmetry and boundary conditions (11 ip with a = T. With the notation of the preceding 
section we can identify b = a = L, i.e. both the position of the local maximim and the 
range of the potential are given by the box length parameter L. Dropping the index 
T we can write the transmission resonance wave function inside the square well as a 
superposition of plane waves, 

ipix) = Fe'^^ + Ge"^^^ , \x\ < L (23) 

where q = ^j2m{E — Vo)/h is real. The symmetry condition |-?/^(— x)p = |?/'(x)p leads 
to FG* = GF* so that 

|^(a;) |2 = + + 2FG* cos(2gx) . (24) 

The second condition in (ITT]) reads ip'{L) = ikip{L). Inserting Eq. ( 123|) we obtain 

G=i— ^e^'^^F. (25) 
q + k 

The condition FG* = GF* then implies exp(4iq'L) = 1 which leads to the resonance 
condition 2qL = rni with integer n. Thus we arrive at the celebrated formula 

^" = ^° + i;^ = ^» + 8;^«- 
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for the transmission resonance energies of a finite square well potential. The number n 
must be sufficiently large to make positive. Inserting Eq. into leads to 

l + 2(-l)"^cos(2gx) ' ("^-^^^ 



q + k 



q + k 



cf. left panel of figure [T]) and in particular 

2 



1 + 



q — k 

q + k 



Q + k^ 



Inserting \il){L)\'^ and the integral 

t 

Jo 



(x)\Mx = L\F\ 



L\F 



[q + kj 

2{q + kf + {q-kf 



2L\F? 



[q + kf 
q^ + k"^ 



(q + ky 

into the Siegert formula (|T8l) we obtain the decay coefficient 

' ^2 



2h 2En q' 



or 



L V q^ + k"^ 



(27) 



(28) 



(29) 



(30) 



(31) 



L S m 2En-Vo 

which is the well known textbook result for the decay coefficient (or resonance width) 
of a finite square well potential (see e.g. [I3]), that is usually obtained by expanding the 
transmission coefficient in a Taylor series around E = En- 



3.2. Delta Shell potential 

As an analytically solvable example for complex resonances in asymmetric potentials we 
consider the one-dimensional delta-shell potential 



V{x) 



+00 



X < 



{h'^/m)X5{x-L) x>0 



with A, L > 



(32) 



consisting of an infinitely high potential barrier at x = and a delta barrier at 
X = L which is illustrated in the middle panel of figure [TJ To find the solution of 
the corresponding Schrodinger equation 



2m dx"^ 
we make the ansatz 



+ V{x) 



^(x) = (E -iT/2)iP{x) 



ip{x) 



I sm(kx) 0<x <L , ^ ,^ n^P 
4kx r andE-ir/2- 



Ce" 



X > L 



2m 



(33) 



(34) 
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for the wavefunction which satisfies the required outgoing wave (Siegert) boundary 
condition for x — ?■ oo. The matching conditions for the wavefunction at a; = L read 

ij{L - e) = ij{L + e) , i:\L - e) + 2X^{L) = ij\L + e) (35) 

where the discontinuity in the derivative is caused by the delta function potential [T3] . 
This leads to 

/ sm{kL) = C , kl cos{kL) + 2X1 sm{kL) = ikC e'^^ (36) 

or 

k cos{kL) + (2A - ik) sm{kL) = . (37) 

The real and imaginary part of the eigenenergy can now be found by numerically solving 
the transcendental equation (I57|) in the complex plane. In the following we show how the 
Siegert approximation method presented in section [2] can be used to obtain a convenient 
approximation in an analytically closed form. To obtain the approximate real part of 
the energy and approximate wavefunction as required in step ([I]) of the method we make 
the following approximation: Imagine that the delta potential at x = L is infinitely 
strong, i.e. the limit A — )■ oo. Then the system is a closed box of length L and the 
wavenumber satisfies kL = nir with an integer n. For a strong but still finite delta 
potential, i.e. \k\ << A, we therefore assume kL = nir + 6L with 6L < and \6L\ ^ 1. 
Inserting this ansatz into the real part of Eq. (1371) and expanding it up to second order 
in SL yields 



^ _ 2AL + 1 
niiL 



Inserting k = rni/L + 5 into equation (13^ yields approximations for the resonance 
wavefunction and the real part of the eigenenergy. An approximation of the imaginary 
part is obtained by inserting these results into the Siegert formula f lT9|) . 



2m/(f |^(x)|2dx m\L J ravr + (5L - sin(2(5L) ' 

where we have identified 6+ = a+ = L and = a- = 0. For a potential with 
A = 10, L = 1 and scaled units h = m = 1 the Siegert approximation yields 
£approx = 4.481 — i0.062 for the most stable resonance (cf. middle panel of figure [1]) 
which is in good agreement with the numerically exact result Sexact = 4.487 — i0.061. 
A treatment of the equivalent problem within the context of the nonlinear Schrodinger 
equation can be found in [8]. 



3.3. Double barrier 

As an example of a numerical problem we consider the double barrier potential 

V{x) = ^x^exp{-Xx^) (40) 
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with Vo = 1 and A = 0.1 so that the position of the potential maxima is given by 
h = 1/v^ ^ 3.16 using units where h = m = 1. In order to calculate the approximate 
resonance wavefunction and real part of the energy as required in step ([ij) of the method 
we solve the boundary value problem given by equation f llip by means of a shooting 
procedure. We choose the cutoff parameter for the potential a = 20 to ensure that 
V{x) ~ for \x\ > a so that the wavefunction is well approximated by a plane wave 
in that region. Starting with initial conditions given in equation ( II ip we integrate the 
Schrodinger equation from x = —a to x = using a standard Runge Kutta solver. 
By means of a bisection method the real energy E is adapted such that the boundary 
condition at x = is satisfied. The decay rate is again obtained by inserting this value 
of E and the corresponding wavefunction into the Siegert formula (1181) . More details on 



the actual numerical implementation of the method can be found in Appendix A 



Table 1. Complex eigenenergies for the three most stable resonances of the potential 
(PO)) calculated using complex scaling (CS) and the Siegert approximation (S). 



n £s 




£cs 






1 0.4601 - 


19.62 X lO-'^ 


0.4601 - 


19.6204 X 


10-^ 


2 1.2804- 


il.70 X IQ-^ 


1.2804- 


i 1.6737 X 


10-3 


3 1.88 - 


i7 X 10-2 


1.8531 - 


16.7240 X 


10-2 



The right panel of figure [T] shows the potential (1400 and the square {ipl"^ of the 
most stable resonance wavefunction calculated with the Siegert approximation method. 
Table [1] compares our results for the three most stable resonances with numerically exact 
results calculated with the complex scaling method described in [T3]. We see that our 
simple approximation yields very good results for the ground state since its decay rate 
is small. The values for the first and second excited states demonstrate that the Siegert 
approximation becomes less accurate with increasing decay rates. A generalization of 
the same problem to the nonlinear Schrodinger equation is straightforward and can be 
found in [7]. 

4. Summary and conclusion 

In this article the Siegert approximation method for calculating complex resonance states 
was presented in an intuitive and straightforward manner which at the same time clearly 
points out the similarities and differences between resonances of transmission coefficients 
(or similar quantities) and resonances in the complex plane (Siegert rsonances) as 
well as the role of the continuity equation in this context. It was illustrated by 
two analytically solvable example problems and a numerical application. The author 
hopes that the present article, in addition to drawing attention to a useful and easily 
applicable computational tool, offers an alternative, rather intuitive point of view for a 
better understanding of complex resonances which complements other, more technical 
treatments. 
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Appendix A. Numerical calculation 

The following commented MATLAB code implements the numerical algorithm for 
calculating resonances of the one-dimensional symmetric finite range potential described 
in section [3731 For paedagogical reasons the program makes use of several global variables 
which give the code a simple structure but make it less flexible and elegant. For 
the same reasons and for achieving compatbility with both MATLAB and the Open 
Source software OCTAVE the integration of the Schrodinger equation is performed by 
means of a straightforward implementation of the classical Runge Kutta method which 
requires a rather high number of grid points. Thus the present code can be made a lot 
more efficient by using more sophisticated integrators like, e.g., MATLAB's 0DE45 or 
octave's Isode. The program can be straighforwardly adapted for other symmetric 
potentials by changing the function V{x) and providing the corresponding position h 
of the potential maximum that confines the wavefunction as well as a suitable cutoff 
parameter a. For a potential open on one side only the boundary condition criterion 
crit must be modified according to equation (12U|1 . 

function [E_res , Gamma, psi , Dpsi , x]= Num_Siegert 

% Calculates complex resonances of a double barrier potential 
% with the Siegert approximation method 
global VO lambda m hbar a Nx 

% Parameters for the double barrier potential: 
V0 = 1; lambda = 0.1; 

% hbar and mass: 
hbar = l; m=l; 

% Shooting method for calculating the resonance wavefunction 
% and real part of the energy: 

a=— 20; % left boundary of integration interval [a, 0] 
Nx=10000; % number of lattice points / steps for integration 

Emin = 0.3; % lower boundary for the energy 
Emax = 0.6; % upper boundary for the energy 



Calculating resonance positions and widths using the Siegert approximation method 12 

% use bisection method to find the real part of the resonance 
% energy of the solution of the Schoedinger equation in the 
% energy interval [Emin, Emax] : 

E_res=bisect ( @ solve _SE , Emin , Emax, 10 " —5) 

% Calculate wavefunction 'psi ' and its derivative 'Dpsi' 
%for the resonance energy 'E_res ' ('x' is a vector of 
% grid points, ' crit ' is (close to) zero in case of 
% resonance (see below ) ) : 

[crit ,psi ,Dpsi,x] = solve_SE(E_res); 

k=sqrt (2*m* E_res ) /hbar ; 
dx=(0-a) /Nx; 

b = l./sqrt (lambda) ; % calculate position of potential maxima 

% Use Siegert formula to calculate the decay rate 
% ( cf . equation (18) in the text) : 

Gamma^hbar " 2*k/m*abs (psi(l))."2./ (sum( abs ( psi ) . " 2 . * ( x>=— b) *dx) ) 
Gamma_haIf=Gamma/ 2 

% plot s olution : 

% normalize wavefunction for plot: 

psi=psi . / sqrt (sum( abs ( psi ) 2 .* (x>=— b ) . * dx) ) ; 

figure ( 1 ) ; 
hold on 

plot (x , abs ( psi ) . " 2) % plot wavefunction 
plot (—X , abs ( psi ) ."2) 

plot (x , V(x) , ' r ' ) % plot potential 
plot(-x,V(x) , 'r ') 
hold off 
xlabel( 'x') 



% 

function [crit ,psi_v ,Dpsi_v,x_v] = solve_SE(E) 

% integrate the Scroedinger equation from x=a to x=0 
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% for a given energy E; Return values: 

% crit : crit=0 if boundary condition at x=0 is satisfied 

% ( cf . equation (11) in the text) 

% psi-V : vector containing the wav efu nction 

% Dpsi-V : vector containing the derivative of the wavefunction 
% X-V : vector containing the grid points of the integration 
global m hbar a Nx 

psi_v = []; % empty vector for the wavefunction 
Dpsi_v = []; % empty vector for the derivative of the 
wavefunction 

x_v = []; % empty vector for the grid points 

k=sqrt (2*m*E) / hbar ; 

% Siegert boundary conditions at x=a: 
x=a ; psi=0.001; Dpsi=i *k* psi ; 

x_v = [x_v X ] ; 
psi_v = [psi_v psi ] ; 
Dpsi_v = [Dpsi_v Dpsi ] ; 

dx=(0— a)/Nx; % Step size for integration 

for j=l:Nx % integrate over Nx steps 
yO = [psi; Dpsi]; 

% Runge Kutta integration step : 

dy_dxO=Schroedinger ( yO , X ,E) ; % call to function 
%that implements the Schroedinger equation 

yA=yO + 0.5*dx*dy_dxO ; 

dy_dxA=Schroe dinger (yA, x + 0.5*dx ,E) ; 
yB=yO+0.5*dx*dy_dxA ; 

dy_dxB=Schroe dinger (yB ,x + 0.5*dx,E) ; 
yC=yO+dx*dy_dxB ; 

dy_dxC=Schroe dinger (yC , x+dx ,E) ; 
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y l=yO+dx * ( dy_dxO +2* (dy_dxA+dy_dxB )+dy_dxC ) / 6 ; 

% Store result of integration step in vectors : 
psi=yl(l); Dpsi=yl(2); 
x=a+ j *dx; % increase x 
x_v = [x_v X ] ; 
psi_v = [psi_v psi ] ; 
Dpsi_v = [Dpsi_v Dpsi ] ; 

end ; 



% calculate criterion for the boundary condition 
% at x=0 ( cf . equation (11) in the text): 
cr i t =( psi *Dpsi ' + psi '*Dpsi) ; 

% 

function dy_dx=Schroe dinger (y ,x ,E) 

% implementation of the Schroedinger equation as a 
% system of two first order differential equations 
global m hbar 

dy_dx = [y(2) ; 

2*m/hbar"2 * (V(x)-E) . * y ( 1 ) ] ; 

% 

function potential = V(x) % Double barrier potential 
global VO lambda 

potential = VO . * . 5 . * x . " 2 . * exp(— lambda . * x . " 2 ) ; 

% 

function [xO] = bisect ( func , xl , x2 , tol ) 

% bisection method for solving a nonlinear equation 

f l=feval ( func , xl ) ; if fl = 0, xO=xl ; return; end 

f 2=feval ( func , x2 ) ; if f2 = 0, x0=x2 ; return; end 
if fl*f2 >= I xl >= x2 

xl , fl 

x2 , f2 

error ( 'No^r oot ^found^due^to ^ initial ^values ^xl , ^x2 

end 



for i=l:100 

if x2 — xl < tol return; end 
xO = 0.5*(x2+xl) ; 
fO=feval (func ,xO) ; 
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if fO*fl <= 

x2=x0 ; f2=f ; 
else 

xl=xO; fl=fO; 
end 
end 

error ( 'No^root ^found ' ) 
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